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Abstract 

Abstract. Eigenvalues and eigenfunctions of Mathieu's equation are found in the short 
wavelength limit using a uniform approximation (method of comparison with a 'known' equa- 
tion having the same classical turning point structure) applied in Fourier space. The uniform 
approximation used here relies upon the fact that by passing into Fourier space the Mathieu 
equation can be mapped onto the simpler problem of a double well potential. The result- 
ing eigenfunctions (Bloch waves), which are uniformly valid for all angles, are then used to 
describe the semiclassical scattering of waves by potentials varying sinusoidally in one direc- 
tion. In such situations, for instance in the diffraction of atoms by gratings made of light, it 
is common to make the Raman-Nath approximation which ignores the motion of the atoms 
inside the grating. When using the eigenfunctions no such approximation is made so that the 
dynamical diffraction regime (long interaction time) can be explored. 

1. Introduction Consider a diffraction experiment in two dimensions, as depicted in Figure [I]. 
A plane wave exp(ifcz) propagates freely in the z (longitudinal) direction and is incident normally 
upon a medium with a refractive index which varies weakly in the x (transverse) direction 

n(x) = no + n\ cos 2Kx (ni^na). (1) 

As the refractive index depends only on x, the wavefunction inside the medium separates, ^(z, x) = 
<fi(z)ip(x), with <p{z) being trivially given by 



4>(z) — exp \ iz y k 2 nl — nj (2) 

with k a separation constant proportional to the transverse energy of the wave inside the medium. 
The transverse behaviour is governed by Mathieu's equation (see Abramowitz and Stegun 1964) 

d 2 tb(x) 

J \ ' + (k + 2k 2 nom cos 2Kx) ip{x) = 0. (3) 

More generally, Mathieu's equation is obtained whenever the 3-dimensional Hclmholtz wave equa- 
tion is separated in elliptical coordinates, useful for, say, scattering from elliptical boundaries. 
This paper is concerned with the short wavelength (semiclassical) regime for which the parameter 
2k 2 noni/ K 2 is very large. One motivation is that the resulting asymptotics are known to describe 
the emergence of interesting classical features such as caustics (singularities of the geometric ray 
theory) which come to dominate the wave field as the wavelength is reduced to zero. The caustic 
structure becomes ever more intricate as the (longitudinal) thickness of the medium is increased 
(Berry and O'Dell 1999). 



Diffraction by a sinusoidal grating has been studied in the context of the diffraction of light by 
ultrasound since at least 1921 (Brillouin 1921) (see Berry (1966) for a review to 1966). More recent 
interest has arisen through the realisation of the diffraction of beams of atoms by beams of light 
(Adams et al 1994). Then n(x) = — V(x)/E, where E is the energy of the atoms and V(x) is 
the potential energy due to their interaction with a standing wave of light (Cohen- Tannoudji et al 
1992, Kazantsev et al 1991) 

v{x) = -^ A2 + rV os2 Kx = ~ v ° cos2 Kx (4) 

with £q the magnitude of the electric field of the counter propagating laser beams which form 
the standing wave, K their wavenumber and A the frequency detuning from resonance. V is the 
spontaneous decay rate for the excited atom and d the atomic dipole moment for the electronic 
transition being used. 

The resonant nature of the atom-light interaction allows for an efficient transfer of transverse 
kinetic energy to the atoms. Combined with their large mass, this means the atoms can attain small 
transverse de Broglie wavelengths and so are rather good candidates to access the semiclassical 
scattering regime when compared with other microscopic particles such as neutrons or electrons. 

The idealised experiment described above assumes that the standing wave laser field has a 
'top-hat' cross-section in the longitudinal direction, switching on at z = 0, and remaining constant 
till switching off at z = Z at which point the atoms propagate undisturbed to a detector in the 
farfield. One might achieve this by telescopically expanding the normally gaussian laser profile and 
physically masking the entry and exit edges to make them sharp (diffraction limited). In this way 
the entry and exit into the laser can be sudden from the point of view of the dynamics of the centre 
of mass of the atom, but still adiabatic from the point of view of the rapid Rabi oscillations of the 
internal electronic states that generate the potential given by Equation (g), see O'Dell (1999) for 
more details. In any case, features such as caustics will still be qualitatively correctly described by 
the simple model given here even if experimental conditions differ considerably. This is because 
caustics are stable to perturbations as guaranteed by the optical catastrophe theory (Berry 1980). 

2. Dynamical diffraction and the Raman-Nath equations 

The time independent Schrodinger equation governing the passage of atoms of energy E — 
h 2 k 2 /2m through an optical standing wave of periodicity itjK is 

a 2 * a 2 * / 2 2mv 2l „ A T n 

^ + ^+^ 2 + — cos 2 (i^)j* = 0. (5) 
The periodic potential suggests an atomic wavefunction of the form 

V{x,z) =e ikz M z ) e2inKx - ( 6 ) 

n— — oo 

An advantage of this decomposition of the wavefunction is that upon exiting the interaction region 
at z = Z, the terms in Equation ([|) represent freely propagating diffracted waves travelling at 
angles arcsin(2ni^/fc) to the z axis with amplitudes A n (z = Z). A detector in the farfield will 
register a diffraction pattern made up of discrete beams with intensities |A„(Z)| 2 . The amplitudes 
satisfy J2n=-oo l^"| 2 = 1- The rest of this paper is devoted to determining the A n (z). 

If the initial kinetic energy of the atoms is thermal then k 3> K and the propagation of the 
atom beam is paraxial. The paraxiality means that the evolution of the A n with z will be much 
slower than exp(ifcz) so when substituting (^|) into (^|) terms containing d 2 A n /dz 2 can be ignored. 
The result is an infinite series of coupled equations 

r)A A 

i^-n 2 A n + -(A n+1 + 2A n + A n _ 1 ) = Q (7) 
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where 



, A S J^. (8) 

The parameter A is equivalent to the parameter 2k 2 noni/K 2 appearing above in Eq. (||) — by 
letting Ti — * 0, and hence A — * oo, one obtains the classical limit. Of course, taking the limit h — > 
is a formal device. In an actual experiment the short wavelength limit is approached by, say, 
making Vo, the interaction between the atoms and the light, as large as possible. This increases 
the depth of the wells of the sinusoidal potential which in turn means there are more quantised 
transverse states. A phase transformation A n — > A n exp(iA£) slightly simplifies the equations to 

" n2An + 2 + An ~ l) = °- (9) 

These are (apart from a straight forward change of variables) the differential difference equations 
introduced by Raman and Nath (1935,1936) to describe the diffraction of light by ultrasound. The 
Raman-Nath (RN) equations are a description of dynamical diffraction, yielding the evolution of 
the amplitudes of the various diffracted beams as the atom wave passes through the light grating. 
In their original paper, Raman and Nath (1935) observed that by ignoring the diagonal term, n 2 A n , 
one obtains simple solutions for A n in terms of Bessel functions. This is equivalent to neglecting 
the transverse kinetic energy of the atoms and the sinusoidal potential then acts only as a pure 
phase grating, see Berry (1966). This is a very successful approximation for short interaction 
times (Sanders 1936, Gould et al 1986, Rasel et al 1995) but dynamical diffraction requires that 
full account be taken of the transverse motion. 

A general property of paraxial systems is that the axial coordinate, in this case the rescaled 
longitudinal distance, Q, plays the role of time. The numerical integration of the RN equations is 
relatively simple since they are first order differential equations in Q with only a single boundary 
condition: A n {C, = 0) = 5 n Q. However, when investigating the behaviour at long interaction 
times it becomes more economic to analyse the problem in terms of the eigenfunctions of the 
scattering potential, which propagate unchanged through the medium. This approach will be 
adopted by seeking eigenfunctions of the RN equations (^) of the form A n = f? n exp(— LE£). 
Corresponding to each eigenvalue E 3 is an eigenfunction consisting of a 'vector' of amplitudes 
(Btoo, ■ ■ • , B 3 _ 2 , B^_ lf Bq,B{, B° 2 , ... , B^), whose elements satisfy 

&Bi = n 2 Bi ± (Bi +1 + BiA . (10) 



This equation defines the tridiagonal RN matrix hamiltonian. For a weak potential (i.e. small A: 
the quantum, non-classical limit) the RN matrix can be approximated by a 3 x 3, or for the special 
case of oblique incidence near a Bragg angle, a 2 x 2 matrix, and analytical solutions for the Bloch 
waves (eigenfunctions) are easy to find (Berry and O'Dell 1998). Here, however, we are interested 
in the opposite limit. 

For the purposes of numerical diagonalisation, a guide to the minimum diffraction order, ±JV, 
at which the RN matrix can be safely truncated for large A is given by 

N = V2A oc 7T 1 . (11) 

This includes only those beams contained within the maximum scattering angle that can be 
achieved classically (Berry 1966). Criterion ( pi] ) becomes exact as A — > oo but at the cost of 
requiring an infinite number of beams. 

The 'physical' derivation of the stationary RN equation ( JlOj ) given here is nothing more than 
a Fourier analysis of Mathieu's equation (||). Equation ([Tof) is the recursion relation satisfied by 
the Fourier coefficients of periodic solutions to Mathieu's equation: the even and odd Mathieu 
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functions (Abramowitz and Stcgun 1964). The eigenvalues E 3 are known as the characteristic 
values. 

3. WKB solution of the Raman-Nath equations 

In the semiclassical limit the dimensions of RN matrix defined by Equation (|l^) are infinite and 
numerical diagonalisation becomes impossible. Dingle and Morgan (1967a, 1967b) (see appendix 
of Berry 1966) and independently Yakovlev (1997), have given a WKB-type solution for the con- 
tinuised RN equations accurate in the large A limit. The idea is to capture the very fast oscillation 
of the Bloch wave (see Figures with the exponential of a slowly varying function. To this end 
one replaces the discrete variable n with the continuous one y 

(12) 

so that when n — > n + 1, then y — > y + (VA) -1 and the discrete amplitudes become continuous 
functions of y, B n — ► B(y). Defining the rescaled eigenvalue 



(3 = 



E 
A 



(13) 



= 0. 



the stationary RN equation (|T(J) becomes 

[fi - y 2 )B{y) + ± [b (y + (VX)- 1 ) + B (y - (VX)- 1 ) 
Following Berry (1966), let 

B(y) = e iSiv} 

where the suggestively named S(y) is analogous to an action. Taylor expanding S(y + (-s/Ap 1 ) 
and S(y — (VA) -1 ) gives a differential equation of infinite order 



(14) 



(15) 



cos S l + 
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( y 2_ ^-1(5^+5-/24+-) 



(16) 



where S m = (VA) 711 d m S/dy m are assumed to be small quantities of order (yA) m . Solving for 
the first derivative one has 



1 OS 



arccos 



VXdy 

Expanding the right hand side (rhs) gives 



(y 2 - ^e" 1 



■S* i /2+S iv /24+-) 



~6~ 



(17) 



y 



s 



1 dS r 2 a, 
—=—— — arccos \y — p\+ 1 

VX dy v/1 - (?/ - 2 

y 2 -P 



S" 



(18) 



(i - (y 2 - P) 2 f 2 8 

which can be solved for S 1 by iteration, yielding to second order 



(19) 



The second term on the rhs can be integrated immediately so that the equation for S can be 
written 



S S3 VA f arccos [y 2 - (3\ dy - i In (l - (y 2 - (3 f 



VI So(y,0)- '-In (1 -(y 2 -(3) 2 ) 



(20) 
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Thus the continuised eigenfunctions take the form (Berry 1966) 

piVA / arccos [y 2 -(3] dy p i\/A So(y,P) 

B(y) = c iS » — = — (21) 

(1-^-/3)2)1/4 (1 _ (y 2_ ^2)1/4 

which resembles a WKB expression. In particular the denominator causes divergences at the 
turning-points 



y = ±^]3±l. (22) 

'Bound' solutions of Mathieu's equation have energies lying between the top and bottom of the 
sinusoidal wells (0 < E < Vq) which translates into the bound eigenvalues occupying the range 
— 1 < f3 < 1. This means that except for the situation when (3 — 1, real turning points for bound 
states are located at 



y± = ±y/fJ + l. (23) 

'Free' states have energies above the wells. When E is sufficiently greater than Vq it is easy to find 
WKB solutions of Mathieu's equation directly in coordinate space (since they have no classical 
turning points and hence no divergences) rather than the momentum (Fourier) space used here. 
Similarly, it is also simple to find eigenfunctions in coordinate space when j3 « — 1, since states 
near the very bottom of the wells are the most localised and see an essentially harmonic potential 
yielding hermite polynomials as solutions. The most interesting situation is for [3 lying close to 
+1. These are the states affected most by tunnelling between the wells and will be examined in 
sections 10-19. 



The remaining 'action' integral, giving the phase of the WKB solution (21), can be calculated 
using the positive turning point y + = y/jT+T as the lower limit (i.e. the zero or reference point of 
the phase) 

fV 2 

So(y+,V,P)= / arccos [y' - (3} dy' 

JVP+I (24) 



y arccos [y 2 - /3] - 2y/l + /3 E Q 



arccos [y — (3] 



1 + /3 



where E(0|m) = yl — m sin 2 9 d9 is the incomplete elliptic integral of the second kind, see 
Gradshteyn and Ryzhik (1965). This expression for the phase is valid for < y < y/(3 + 1. For 
perpendicular incidence the phase is symmetrical about y = 0, and so only this half-range is 
required. For values of y greater than ^/ (3 + 1 the phase is purely imaginary, but with care ( |2~i| ) 
still gives the correct answer. 

4. Single and double wells in momentum space 

Whilst it is useful to think of Equation (|2|) as a WKB type expression it does have some 
unusual features due to its unorthodox derivation from a difference equation. Usually, for the 
Schrodinger equation 

^> + ^(,) = (25) 

where p(q) is the momentum, one has the approximate WKB solution (Berry and Mount 1972), 
valid for small h as long as one is not too close to the turning-points {p(q) =0), of 



^wkb = ^r= exp p(q')dq'J (26) 

where +/— refers to right/left travelling waves. Equation ( pl| ) is actually for the momentum 
space wavefunction, but to keep the analogy with the familiar coordinate space WKB solution (|26| ) 
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simple, Equation ( pl| ) will temporarily be treated as though it is a coordinate space expression. 
Thus, terms such as 'momentum' will refer to functions playing the analogous role to p(q) above. 
In particular, what is peculiar about solution is that the 'momentum' function appearing in 
the amplitude and phase are different. The two momenta, 



Pl (y,0) = ^Jl-(y 2 -0) 2 (27) 

and 

p 2 (y,/3) = arccos [y 2 - /?] (28) 

coincide for — > — 1, but are quite different when — > 1. Examining Figure || one notes that 
the momentum appearing in the phase, p2, is that exhibited by a particle in a simple well. The 
amplitude momentum, pi, however, corresponds to a particle in a double well although, except 
for values of greater than one, the particle has enough energy to move between the two wells. 
It is the classical turning-point structure that is of paramount importance (see Berry and Mount 
(1972) for a review of the WKB procedure), so the different local values of the two momenta give 
similar behaviour when < 1. However, as the turning-point structure of p\ changes from two to 
four, at = 1, one can expect a qualitatively different response. 

This retrospective observation, that the RN equations in some sense describe a wave in a 
single/double well, will be exploited to find uniform solutions in section 7. 

5. The Bohr-Sommerfeld condition 

For states in a well, single-valuedness of the wavefunction dictates that only certain discrete 
energies are allowed. These eigenvalues/characteristic numbers ensure the integral of the WKB 
phase, Equation (p4|), from one turning-point to the other (that is, the integral across the clas- 
sically accessible part of the potential well), correctly matches the oscillating part of the WKB 
wavefunction onto (asymptotically) exponentially decaying parts of the wavefunction (that tunnel 
into the classically forbidden sides of the well). 

The action right across the well, given by setting the upper integration limit of Equation (|2~i| ) 
equal to —\J0 + 1, is also equal to, for perpendicular incidence, twice the value found by integrating 
only halfway, to the midpoint of the well at y = 0, 



2S (y+, 0, 0) = - + E Q arccos [- 



■0\ 



1 + 



(29) 



The Bohr-Sommerfeld condition then states 



2VA \S (y+, 0,^')| = (j + ^Tr, J = 0,1,2,3... . (30) 

The root of this equation (which must be found numerically) for each value of j gives the 0° 
eigenvalues. This expression also gives the number of bound states that exist once a value for A 
has been chosen. The term 'bound' here refers to the coordinate space states that energetically lie 
below the maxima of the sinusoidal potential. In momentum space all the states are trapped in 
a well. The most energetic bound state, labelled j ma x, has the value of which is closest to one. 
When A is large the eigenvalues lie very close to each other, and in particular, 

lim j ™* = 1 (31) 

A-+oo 



and since E 1) = 1, then 



hm J m ax = 7T- ( 32 ) 

A— >oo 7T 2 
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Although this gives the total number of possible bound states, not all of them are necessarily used 
in the superposition which gives the diffracted wavefunction. The superposition coefficients, which 
establish the relative contribution of each Bloch wave to the total wavefunction, are determined 
from the overlap of the initial plane wave (A n (C, = 0) = 6 n o) with the Bloch waves. The overlap 
integral is trivial in momentum space. Each coefficient is given by the value of the corresponding 
Bloch wave at y = n = 0. Thus only the even Bloch waves are excited. 
6. Real eigenvectors and normalisation 

Since the stationary RN equation (|l4| ) is real it is always possible to find real solutions. This 
is achieved by constructing a superposition of the two independent solutions: the right and left 
travelling waves of Equation (^6|) , which correctly matches the exponential decay into the classically 
forbidden regions (Berry and Mount 1972). In the classically allowed region the single well has the 
solution 

Bwkb(!/, 0) = /4 cos (VA S (y+, y, (3) + j) (33) 

(l - (?/ - P?) 

where So is given by Equation (|24|), and ftf(0) is a normalisation factor 



Normalisation of the discrete amplitudes B n requires that Y^=-oo l-^«| 2 = 1- When moving 
from a summation to the integration in the continuous variable y, care must be taken to include a 
factor of \/A which comes from the definition ( |l2| ) of y, so that 

oo 

E - 



dn — > VA / dy. (34) 



And so normalising the eigenvectors requires the evaluation 



of 



2 



BWb(v) dy = 1- (35) 

When A becomes large the exponential decay of the wavefunction into the sides of the well becomes 
very rapid and one can ignore these contributions, so the integral is taken to be just that between 
the two turning-points of the classical motion. Further, the oscillation of the wavefunction is 
assumed to be very rapid (actually the normalisation factor derived in this way works well even 
for the ground state Bloch wave which has the shape of a Gaussian, i.e. is non-oscillatory) in 
comparison to the slow variation of the square of the amplitude. And so, without much loss of 
accuracy, the cos 2 term can be replaced by its average value of one half. Although the amplitude 
diverges at the turning-points this divergence is still integrable. Thus, one takes 

VA /■»+ dy „ J_ 

2 J y _ Vi - (y 2 - P) 2 ~ ^ ( j 

which gives 



TV 2 w ^ (37) 



Vak(im) 



where K(m) = (1 — msin 2 ff)^ 1 ! 2 dO is the complete elliptic integral of the first kind (Grad- 
shteyn and Ryzhik 1965). 

7. A Uniform approximation for the Raman-Nath equation 

If the total diffracted wavefunction was calculated as a sum over the WKB Bloch waves then the 
diffraction pattern would contain spurious divergences at the turning points of each eigenfunction. 
This problem can be overcome using uniform approximations which give smooth and uniformly 
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accurate eigenfunctions. Uniform approximations are based upon the idea that one can express the 
solutions to an unstudied differential equation in terms of those of a well known, studied, differential 
equation provided the two share a similar transition point (classical turning point) structure. A 
review of the uniform method can be found in Berry and Mount (1972). Since the method is central 
to the following calculations, it is reviewed in Appendix 1. Somewhat non-standard, and to the 
best of the author's knowledge, novel uniform approximations for the momentum space Mathieu 
functions will be described in this and the following sections. 

If one were to treat Mathieu's equation directly in coordinate space then the infinite number 
of turning points due to the periodic well structure of the potential make a semiclassical analysis 
more complicated, see Berry (1971). Approximations local to one or two turning points (e.g. the 
expansions due to Sips, see Abramowitz and Stegun 1964) must be carefully joined together to 
give the complete solution. By using an imaginary coordinate the Mathieu equation is converted 
to the modified Mathieu equation, Abramowitz and Stegun (1964), with a hyperbolic cosine po- 
tential, and so a uniform approximation for a single well could be applied (see for example Ancey, 
Folacci and Gabriclli (2000), which draws on the seminal contribution by Olver (1954), for uniform 
approximations to the modified Mathieu equation in coordinate space). However, for the type of 
diffraction problem treated here a uniform approximation made directly in Fourier space is more 
useful for two reasons. Firstly, if it is the farfield diffraction pattern that is required, then this is 
the Fourier space description and no further transformation is required. Secondly, for the incident 
plane wave initial condition considered here, the calculation of the superposition coefficient for 
each eigenfunction contributing to the total wavefunction is, as described earlier, trivial. 



When constructing a uniform approximation to the continuised stationary RN equation (14) it 
is not obvious what 'momentum' (the pretense of being in coordinate space rather than momentum 
space will continue to be maintained) function (i.e. p(q) = \J E — V(q) in Equation fl26|)) to use. 
The WKB-type solutions ( plj ) continue to play an important role since they suggest using the 
'momentum' functions p\ and p 2 of Equations (|27|) and (|2"8|). It turns out that, as in the WKB case, 
matches to the actual momentum space Mathieu eigenfunctions are found when pi is used for the 
'amplitude' part of the uniform approximation and pi for the 'phase' part. The mapping function 
should be based upon p%, which has the structure of a simple well. The simplest comparison 
equation for a well is 



da 2 



+ (t-o* 



= 0. 



(38) 



The parameter t depends on the energy. The 'equivalent points' needed for Equation ( |107|) are 
chosen to be the turning-points. This of course im medi ately satisfies the requirement that the 
zeros of T and \ correspond (See Appendix). Using (107), the integral across the well gives t, for 
then 



[ V+ VX P2 (y,f3) dy=-2VXS (y+,0,f3) = [ + ^ y/t - a 2 da 



tn 
T 



(39) 



where So is given by Equation (g4 
function a{y) can be found from 



Once t, which is a function of /3, is known, the mapping 



v7 



V 1 1 — a 2 da 







a 


\ 


arcsin 


LVtJ 
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t 



(40) 



Clearly this step must be executed by numerical root finding for each value of y which is required, 
i.e. those spaced at l/VA intervals which are the angular positions of the diffracted beams. 

There are two standard forms of the parabolic cylinder equation (Abramowitz and Stegun 1964) 



dg< 



±o9 = 0. 



(41) 
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The well Equation ( Bq ) corresponds to taking the upper signs. The simplest independent solutions 
to the parabolic cylinder equations are an even and an odd power series. However, combinations of 
these two power series lead to another two independent solutions, the Whittaker functions, which 
for large g decay or grow exponentially. The Whittaker functions have the correct properties to 
match the exponential tunnelling of the physical solution into the sides of the well. The Whit- 
taker solutions to the well equation (|38| ) are Z3( t _ 1 )/ 2 ( cr V / 2) and Dr t _iy2(— cry/2). The standard 
theory for the potential well uniform approximation would then predict the form of the Bloch 
wavefunction, correct for all y, to be 



B(y) = ^uniform = \n 2^ (JjL 



*(/3)/4 



-1/2 



D 



(t(/3)-l)/2 



{-a(y)V2 



(42) 



and it is noted that for perpendicular incidence the choice of +er or —a makes no difference (only 
even eigenvectors are excited) . The prefactors ensure that this expression has the same asymptotic 
behaviour when a, y — ► oo as the WKB solution (2l|). 

Inserting the Bohr-Sommerfeld condition ([30) into the equation for t, (|39| ) above, gives the 
value of t which corresponds to the j th eigenvalue 

t = 2j + 1. (43) 

When the index of a Whittaker function is an integer, as here, it takes on the more familiar form 

D.j(aV2) = 2-^ 2 H 3 (a)e-' j2/2 (44) 

where Hj is a Hermite polynomial. However, for the higher Bloch waves (e.g. j greater than 20) it 
is more convenient to use, for reasons of speed of computation, the Airy function approximation 
to the Whittaker functions (see Abromowitz and Stegun 1964). 
8. Modifying the amplitude 

The application of the uniform method to the WKB approximation to the RN equation as given 
above requires some adjustment. The existence of two momentum functions means the amplitude 
term of Equation (E||) , 



-1/2 



t-a 2 



arccos 2 [y 2 — (3] 



1/4 



(45) 



does not match the WKB behaviour (see Figure ||), since in that expression it is p\ that appears 
in the amplitude. This disparity is removed by the substitution 



t-a 2 



arccos 2 [y 2 — (3] 



1/4 



t-a 2 



1/4 



1 - (y 2 - PY 



(46) 



The uniform approximation is formulated so that T(a) and x{v) approach zero together so that 
the divergence inherent in the WKB solution is tamed. For the adjusted amplitude to give sensible 
answers this swapping of the momentum expressions must still lead to the correct behaviour at 
the turning-points. For j3 < 1, both p\ and P2 have the same zeros, and crucially for the uniform 
approximation they go to zero in the same way, namely as the square root of the distance from 
the zero. 

L'Hopital's rule can be used to find the limiting value of the amplitude (|4^) at the turning-point. 
After some calculation one finds 



lim 



t-a 2 



1 - (y 2 - PY 



lim 



lim 



sft 



2VTTP 



2/3 



(47) 



9. Comparison with the purely numerical calculation 

Some pictures will now be used to compare the uniform method with the results of numerical 
diagonalisation (which can be taken as the 'exact' result). Figures |]-f7| show a selection of Bloch 
waves with the uniform calculation shown as a solid line, though as before only the discrete values 
of y corresponding to the diffracted beams were used. The dots are the numerical data. 

Figures ||-[l2] give the square of the modulus of the total wavefunction (i.e. the diffraction 
pattern), found by summing the Bloch waves, for a selection of depths (i.e. thicknesses of interaction 
region). As mentioned at the end of section 5, the superposition coefficients in this sum are given 
by the value of the particular Bloch wave at y = 0. Only the bound states were used in both 
the numerical and uniform calculations, an approximation which is increasingly accurate in the 



classical limit. The intensity shown in Figure 



clearly displays (the square of) an Airy function, 



which is well known to be the wavefunction associated with a fold caustic (Berry 1981). Caustics 
are the foci of the diffracted wavefield — classically they diverge and hence come to dominate the 
diffraction pattern in this limit. Figures are the quantum equivalent of vertical slices through 
the classical ray trajectory picture, Figure lb, of Berry (1999) (the depth, £, is the same classical 
unit denoted by x in that paper). The caustics are the envelopes of families of rays (each ray 
corresponds to a classical atom) which oscillate back and forth in the potential. 

If the potential was harmonic then all the atoms would be focussed at y = when Q = mix 
(to = 0, 1, 2 ... ) and this behaviour can be seen in Figure^ for to = 1. However, the anharmonicity 
of the sinusoidal potential means that the foci are imperfect and are smeared out into cusps with 
fold caustic arms. With increasing depth, successive oscillations of the atoms each introduce a 
new cusp which develops into a fold as it moves outwards through the pattern with increasing 
The result is an increasingly complicated interference pattern between sucessive Airy functions as 
shown in Figures [Toj— [l~^, but usually with a few fringes of an outlying Airy function visible. The 
proliferation of caustics for increasing £ is described only by dynamical diffraction and lies beyond 
the phase grating (or so-called 'Raman-Nath') approximation. 

10. The problem of the separatrix 

As has been already been noted, when the eigenvalue (3 approaches 1 the solution to the RN 
equation changes its nature. Classically, the motion of a particle in the sinusoidal potential goes 
from being trapped in a single well (libration) to being free (rotation) when its energy passes 
through the separatrix at /? = 1 from below. At the same time the number of turning points of the 
'momentum' function pi jumps from two to four and in so doing p\ becomes qualitatively different 
from p2 — as apparant from Figure || when the central dip breaks through the zero line. These new 
turning-points occur at 



V 



±Vf3^1- (48) 



For small values of A the actual divergences due to these turning-points can fall between the 
diffracted orders and go unnoticed. As A is increased this is no longer the case and the divergences 
become clearly defined as the classical distribution emerges. It is emphasised that these divergences 
only affect those eigenvectors with eigenvalues close to the separatrix. Careful examination of the 
picture of the 200 th eigenstate for A = 12500, reveals the first hint that the uniform approximation 
has a defect when (3 begins to approach 1. The remainder of this paper is concerned with these 
eigenstates close to the separatrix. Although only forming a small fraction of the total eigenstate 
sum giving the diffracted wavefunction, they are perhaps the most interesting states as they contain 
the very fine corrections due to tunnelling between the coordinate space potential wells in the 
semiclassical limit. 

The uniform approximation used so far was not designed to handle the new turning-points. 
The momentum function pi used in the mapping relation ( |40| ) contains no information concerning 
the new turning-points. The amplitude and phase functions no longer act in concert. One way 
to proceed is by a transformation upon the RN equation which results in both the momentum 
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functions have turning-points at y = V '(3 — 1. Denning 



B n = (-1)"C„ (49) 
the stationary R-N equation (|lj) becomes 

1 



{y l - 0)C n + -(C n+1 + C n _i) = 0. (50) 



from which one obtains an altered action 



^ = VA arccos [/? - y 2 } + i - (/ l/\ 2 (51) 
leading to the WKB formula 

„±i%/A / arccos [0-y 2 ] dy „±i\/A S (y,f3) 
C(V) = T7J~ = — 777 • (52) 

(l-(y2_ ^2)1/4 (1 _ (y2 _ ^2)1/4 

The amplitude is the same as before, giving the two sets of turning-points, but the phase momentum 

p 2 (y, (3) = arccos [13 - y 2 ] (53) 



now has its turning-points at y = ±\/ (3 — 1 as promised. The cost is the loss of the turning-points 
at y = ±Vp + 1. The effect of the transformation is in swapping the roles of inner and outer 
turning-points. Figure^ shows that whereas p 2 has the momentum profile for a well, p2 has that 
of barrier which the particle has enough energy to surmount. 

Since the phase momentum functions p 2 and p2 only describe one set of turning-points each, one 
is forced into employing two separate transitional uniform approximations for each eigenfunction 
when one is close to the separatrix. One transitional approximation covers the inner turning points 
and the other the outer. Both are valid in the intermediate region where they smoothly join. 

11. The parabolic barrier equation 

The inner turning-points require a transitional approximation for a (smooth) potential barrier. 
A suitable comparison equation is 

g + (^ 2 )0 = O. (54) 

For t > 0, the Bloch states are more energetic than the central potential barrier, and classically 
one has transmission above the barrier. This is referred to as the underdense case. Making the 
change of variables 

t = it (55) 
a = -?=e^ 4 (56) 



one is lead back to the equation 



d z (j) ft 



da 2 V 2 



= (57) 



which is the same as the first (upper sign) parabolic cylinder equation ( pd| ) when the identifications 
g — a and a — i/2 are made. 

The appropriate solutions for the barrier top are not the Whitakker functions since a parabolic 
barrier does not give an exponentially decaying wavefunction for large g. Instead, for a barrier, 
the basic even and odd power series solutions, which will be referred to as 0i(a,g) and 02(a,<?) 
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respectively, are the correct choice. Close to the barrier top a is small and the power series solutions 
are most conveniently expressed in terms of the confluent hypergeometric functions 



= e" 9 /4 1F1 



@i(a,g) 
83(0,3) = 3 e-* 2 / 4 



a 1.1. 5^ 

2 + 4' 2' T 

a 3 3 a 2 



(58) 
(59) 



As discussed previously, the boundary conditions mean that only even eigenfunctions are of interest 
here. Thus, the underdense inner turning point transitional approximation will be based upon the 
even power series 



9i 



17Z J 4 



(60) 



12. The action for an underdense barrier 

The underdense barrier does not induce any real turning-points (though of course the proximity 
of the turning-points to the real axis gives the deviation of the WKB amplitude from the true value) 
so the natural choice of reference point from which to integrate the phase is y — a = 0. One finds 



5 (0,y,/3) 



P-V' dy' 



= y arccos [/3 - y 2 ] + 2\\J\ - f3 E ( ~ arccos [/3 - y 2 ] 



1-P 



(61) 



E ( - arccos [0\ 



1-/3 



Although it appears that this action contains an imaginary piece this is actually not the case. 
Strictly, the well known transformations (see Abromowitz and Stegun 1964) should be applied to 
the elliptic functions so that their parameters lie between zero and one (the parameter used above 
tends to infinity as f3 — > 1). When this is done the action 5o is explicitly real. However, the 
transformations produce more complicated expressions so will not be applied here. 

To find the value of t, which was previously given by the integral across the well, one must now 
integrate up the imaginary axis between the points 



y ± = ± v / ^T=±iv / l 3 ^- 
The equivalent points for the underdense barrier comparison equation ( |54| ) are 

a = ±h/t. 

Letting y = iv and a — k, t is implicitly given by 



2iVX 



arccos 



[(3 + v 2 ] dv = 2i Vt 2 
Jo 



q 2 cfc 



which, in a similar to fashion to before, results in the condition 

2 



4v / A v / W? E - arccos 



1-/3 



The mapping function between a and y implicitly giving a{y) is 
VAS o {O,y,0)= I y/t + a 2 da-- 







^arccosh 









a 



(62) 
(63) 

(64) 
(65) 

(66) 
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which, together with the value of t, gives the transitional approximation to the wavefunction 



B(y) = (-l) n C(y) = (" Intranational 

cc (-1) 



t + CT 2 



1/4 



72 



1 1 



4' T 



(67) 



This is a real function for real t and a, which is valid from y = and almost all the way to the 
outer turning-point, breaking down close to it because it is only set up to deal with the inner 
turning-point. The constant of proportionality will now be obtained by matching the asymptotic 
behaviour of this function to the WKB solution somewhere between the two transition points. 
13. The asymptotics of the barrier transitional approximation 
The confluent hypergeometric function has well known asymptotics. When \<j\ is large 



li 7 ! 



t 1 1 



4 4' 2' 



r(§) 



_ e -M-it+l)/4f_ ia 2U«-l)/4 ( 1 



r(|) 



r(|-i|) 

where Y is the Gamma (factorial) function. So 
t 



(it 2 -At- 3i) 
+ 16^2 

(-it 2 — 4i + 3i) 



(9 



^ (68) 



9i 



-1-, V2<re" 



i7r/4 



1*1 



,t 11 . , 

-1 — I — ; — ; —icr' 
4 4' 2' 



r ( i ) (-- 2 ) 



/_. 2V'/ 4 

-7T(t+i)/4 p i^/2 l~ 10 " J r ln[l+(it 2 -4t-3i)/(16o- 2 )] 



72. 



-i. 2 ) 



-it/4 



r(i-4) 



^ln[l-(it 2 +4i-3i)/(16<T 2 )] 



(69) 



which conveniently reduces to 



|r(| + i|) 



.^-1/2^^/8-7(4^) cos f — - In cr 7 — — — Ar 







7T 


Hi 










~ 8 



(70) 



When A is large enough cr quickly takes on large values for even modest sizes of y, and so the 
confluent hypergeometric function attains its asymptotic form in the region between the inner and 
outer turning-points. It may then be compared to the WKB solution ( |52| ) for the transformed RN 
equation. The left and right travelling WKB waves (p2) are combined to give a real solution 



-Bbarrier(y) — 



(l - (y 2 - (3f) 



1/4 



(VAS (0,y,/?)7M/3)7^7n/ 



(71) 



where the (—1)™ factor has been incorporated into the phase of the cosine as y/Airy. To enable a 
direct comparison, the phase of the cosine of Equation (|70]) should also be augmented by the same 
quantity. The real phase angle /i(/3) for this parabolic barrier approximation (which for a simple 
first order turning-point, due to a linear potential, is equal to tt/4) will this time be determined by 
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consistency with the asymptotic solution (fTCf). In order for a comparison to be made, the action 
<So(0, y, (3) appearing in the WKB solution must be written in terms of (a, t), which is accomplished 
through Equation (pq). Expanding the rhs of (p6[) for a ^S> t, one has 





/ o- 2 


^arccosh 






V 1 + T 




t t t a z t t z . 

-W--lnt + -k2 + T + -+0 L] (72) 



implying that 



t t t 
fj, = - hi t In 2 Arg 









^)] 







(73) 



Figure |lj demonstrates that this expression for is correct by comparing the WKB solution ( |7l| ) 
containing it, with the fully numerical calculation. The value of A is reasonably small so the WKB 
solution diverges only very slightly from the correct value. 

The exact solution to the parabolic cylinder equation has thus contributed to the evaluation 
of the phase of the WKB solution. On the other hand, the WKB solution indicates the necessary 
modifications needed for the amplitude of the parabolic cylinder equation so that it becomes the 
correct transitional solution to the particular problem being dealt with. Equating the ampliutdes 
of Equations @ and @), one finds Equation (|67]) can now be updated to read 



barrier 



(y) = (-1)"A-1 



= 7ri/8 



t + (7 A 



1/4 



2T 



i - (y 2 " PY 



72 



1F1 
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4' 2' 



(74) 



The parabolic transitional approximation is compared to the fully numerical result in Figure 
|l5] A). At first sight the match does not seem too good. The reason is that the normalisation uses 
the WKB amplitude factor, which diverges at the turning-points. When there are only the outer 
turning-points this method seems to work (see Figures |^-^|) since the divergences are narrow enough 
to not produce too significant a contribution. However, the appearance of the inner turning-point 
divergences close to the separatrix energy now means the normalisation factor is significantly over 
estimating the magnitude of the wavefunction, and thus reduces the magnitude too much as shown. 
With relatively little effort one can numerically normalise the uniformly calculated eigenvectors by 
summing the discrete amplitudes, and when this is carried out the match, shown in Figure [If] B), is 
exceedingly good. This illustrates that it is only the normalising factor which is at fault. Figure |l5| 
B) further illustrates that the barrier transitional approximation, Equation (fri]), is correct nearly 
throughout the entire momentum range — only breaking down close to the outer turning-point. 

14. Calculation of the eigenvalues close to the separatrix — a modified Bohr- 
Sommerfeld rule 

There is a slight complication to the calculation of the allowed values of (3 close to the sepa- 
ratrix which needs to be highlighted. When comparing the values of [3 obtained by the numerical 
diagonalisation technique with those obtained via Equation (j30|), the two differ when j3 grows very 
close to one. Somehow the derivation of the basic WKB solution ( pl| ) has failed to capture the 
full behaviour of the P2 function — perhaps it should now afterall contain two turning-points, not 
one, and so match the structure of the amplitude p\ term? (Implying the transformation ( |49| ) 
of the phase momentum is more than a device.) From the point of view of the eigenvectors this 
can be overcome by replacing the previous single uniform approximation with two transitional ap- 
proximations when (3 approaches one; the parabolic transitional approximation to cover the inner 
turning-point, and an Airy function approximation for the outer turning-point (since this remains 
a simple first order turning-point). However, to calculate the allowed values of the action which 
corresponds to the bound states, one needs some expression which is valid throughout the entire 
region which joins the two turning-points. 



14 



The general procedure for finding the action across a classically allowed region which separates 
two arbitrary types of turning-point employs two transitional approximations which are each valid 
at one end of the region, but these must be correctly joined. The quantised values of S, and hence 
/?, are those which correctly match the two somewhere in the region of mutual validity. 

The matching is most easily accomplished using the asymptotic forms for the two transitional 
approximations — which are of course their WKB approximations. In the region between the two 
turning-points one thus has 



-jjj cos (VAS o (0, y, 13) + n(fi) + VXiry) 

(l-(2/ 2 -/3) 2 ) V 1 



(l - 0,2 - /J) 2 ) 



(75) 

cos (Vas (a/T+A y, (3) + J) 



which implies that 



VXS o (0, y, 13) + n{(3) + VAvry = VASo(y/l + P, V,P) + \ (76) 



modulo 2-7T. 

The method described above works in conventional situations with WKB expressions developed 
from (continuous) differential equations. Once again however, the approach has to be modified for 
the RN equation — whilst successful for the single well, as soon as the inner turning-points begin to 
approach the real axis even the matching of the two transitional approximations runs into trouble. 
The reason is that the continuous descriptions embodied above by Equation ( f75| ) do not match at 
all. Only when they are evaluated at the discrete points corresponding to diffracted beams do they 
match. The transformation (^) has produced two different equations whose continuised WKB 
expressions only respect their common origin at the discrete level. It is then a surprise to find 
that at the correct (characteristic) values of (3 the discretely evaluated expressions on either side 
of Equation (|7^) are in perfect agreement for all y. Both are identical in each other's supposedly 
exclusive region of validity. This is rather curious, but the characteristic values of f3, which one is 
able to predict by correctly matching the discrete points of the two WKB expressions, demonstrate 
that it is correct. 

Due to the simultaneous validity for all y, the most sensible point to choose to match the two 
solutions is y = 0. The correct matching condition for even eigenstates becomes one of 

cos( A i(/3))-cos(VAS' (j/ + ,0,/3)+7r/4) = (77) 
cos(M(/3))-cos^VA5o(y + ,0,/3) + 57r/4) = (78) 

the choice depending on whether the terminating Airy function has its peak above or below the y 
axis. In fact, successive even eigenstates alternate between the two conditions. When using ( |77| ) 



and (78), it is necessary to express fj,, which is in the first instance a function of t, see Equation 



3h, as a function of (3 through the definition of t (Equation (|65|)). A further subtlety concerning 



the use of (77) and ( |78| ) is that close to each of the characteristic values there is another zero 
which does not correspond to an eigenvalue. The correct zeros are those through which the l.h.s. 
of Equations (|77]) and ([78]) have negative gradients. Table [TJ compares the values of the top eight 
bound eigenvalues for A = 12500 as calculated by the different methods which have been outlined 
so far. Clearly the modified method gives excellent agreement with the true value and is superior 
to the regular Bohr-Sommerfeld scheme when close to the separatrix. The remaining error between 
the modified method and the true value becomes smaller as A — > oo. This is further emphasised 
by the last entry on the table which is the last bound eigenvalue for A = 250000. The Bohr- 
Sommerfeld method predicts only 898 even bound states whereas the modified method accurately 
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finds the value of the 900 th . The subtle behaviour of the eigenvalues near the separatrix can be 
physically attributed to the (semiclassically) exponentially small corrections due to wavefunction 
tunnelling. This phenomena has been recently discussed in a similar context by Waalkens, Wiersig 
and Dullin (1997) and by Sieber (1997). 

15. The Airy transitional approximation 

As has already been pointed out, to obtain the complete wavefunction correct for all y one 
must join the parabolic barrier approximation (|7^ ) to another transitional approximation which 
covers the outer, first orde r, turning-point. The comparison equation is given as an example in the 



Appendix (Equation (|109| )), and choosing the reference point as y = y+ = \/l + (3, the mapping 
function a(y) is given by 

f -1 M 3/2 Hy<VT+P; (a<0) 
AS (y+,y,f3) = { (79) 
[ §kr 3 / 2 ify>yrT?; (a>0) 

since the expression given for So, Equation (|24|), is positive imaginary when y > y + , and negative 
real when y < y + . 

The well known asymptotics of Ai(cr) when a 3> are 

M(a) ~ ±- a" 1 / 4 e"i ff3/2 (80) 
and so the Airy transitional approximation becomes 

transitional = B Ahy (y) = 2tvM ( ^ „ \ Ai ( CT (y)) . (81) 

\l-(y 2 -(3) J 

16. The free eigenstates 

As emphasised previously, 'free' is a description which refers to the (actual) configuration space 
situation of states having transverse energies greater than Vq. In (actual) momentum space there 
are no free states, the classical bounding of the maximum being set by the initial transverse mo- 
mentum plus whatever the atoms can extract from the potential — which depends on the (actual) 
configuration space point, but has a maximum of \/2mVo. Thus, even for (3 > 1, one expects 
caustics in (actual) momentum space. One sees why the free eigenstates are quantised and not 
continuous in energy. Somewhat perversely, the states which are free in (actual) configuration 
space, sit in a double well in (actual) momentum space, and so the central barrier is now over- 
dense — meaning that classical transmission is forbidden. For perpendicular incidence, the free 
'states' are classically inaccessible, so their contribution to the eigensum of states forming the total 
wavefunction is exponentially small. 

For states with (3^1, the problem is most easily solved using the WKB technique in (actual) 
configuration space, since there arc no turning-points to contend with. Constraining the discussion 
to perpendicular incidence means however that only those states with (3 a little greater than one 
need be calculated, so the 'close to the separatrix' treatment of the preceeding sections must be 
generalised to encompass (3 > 1. Since the essentials of the application of the uniform method to 
the Raman-Nath equation have already been conveyed, the following treatment is intended to be 
more of a 'recipe' than a detailed account. 

The overdense barrier equation will be taken as 

+ (- 2 -i)0 = O (82) 

with t a positive quantity. The connection with the parabolic cylinder equation (^) is made with 
the aid of the transformations 

a = -i~ (83) 
g = V^ae™^. (84) 
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To remove any ambiguity regarding the phase momentum function p2 for the barrier, it will be 
written as 



P2 = arccos 



[P-y 2 



i arccosh [(3 - y 2 } if < y < y/JT^T 



(85) 



7r — arccos 



where the central barrier lies between ±^//3 — 1. The actions generated from these momenta, using 

2 



y = \/[3 — 1 as the reference point, are 



q 
and 



( y/0 - 1, y, f3) = i ( y arccosh [(3 - y 2 } + 2iV J) ~ IE ( - arccos [p - y 2 ] 



1 



1-/3 



(86) 



- 1. V, P) = *V + 2^/3+ IE Q arccos [y 2 - /?] 



- 2V^TTe ( \ 



1 + p 



l + p, 

y arccos [y 2 — 0\ 



(87) 



As before, the comparison equation (82) gives rise to the mapping function by setting 







a 




arcsin 


[Vt\ 







and 



5, 







/ \/o 2 -tda = 



(89) 



In particular, the 'barrier integral' which fixes the value of t once (3 is known, can this time be 
conducted along the real axis, and gives, using (|8q) and 



2iVAV/? - IE Q arccos [(3] 



1-/3 



The correct solution to the barrier equation is still the even power series 

e 1 (o,5) = 1 (-i^V^ae^ 4 
and so the transitional approximation for the overdense barrier becomes 

B(y) = (-l) n C(y) = (-Intranational 



(90) 



(91) 
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(92) 



17. Asymptotic matching to the overdense WKB expression 

The transitional wavefunction differs by a few sign changes from the underdense case, and for 
large a these produce the modified oscillatory behaviour: 



|r({- + if) 
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(93) 
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Expanding the rhs of Equation (189) for a t gives 




arccosh 



Vt. 



• In a 
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(94) 



from which one deduces the unknown phase angle /i, appearing in the WKB approximation for the 
overdense barrier (see Equation (f7l|)), to be 



/i = — In f H — In 2 H h Arg 
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(95) 



18. The overdense eigenvalues 

Once again fi can be successfully employed in the accurate determination of the eigenvalues (3. 
Following the empirical observations from the underdense case, the WKB expression emanating 
from the outer turning-point and that from the inner turning-point are matched at a point y 
corresponding to one of the beams. This time the choice of y — is not available since only 
the phase for the WKB approximation outside the barrier is known. The next most obvious 
choice is either the inner or outer turning-point since there the phase of the WKB expressions 
are simplest, but in general these classically determined points will not fall on a diffracted beam. 
Selecting a random beam, y = m/vA, with m an integer, giving of value y lying between the two 
turning-points, will suffice. The condition giving the permitted values of (3 for even eigenstates 
then alternates between 



- cos ( VASo ( VIT/3, + J 



(96) 
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and 
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(97) 



= 0. 



Both of these equations have zeros which do not correspond to the eigenvalues, the correct ones 
being those for which gradient of the lhs' are positive (this is the opposite of the underdense case). 
As before, the accuracy which is achieved gives confidence to the method: for A = 12500 the first 
two free eigenvalues given by numerical diagonalisation are (3 — 1.003356 and (3 — 1.012155, for 
which this WKB matching technique gives (3 = 1.003358 and j3 — 1.012156 respectively. 
19. The overdense eigenvectors 

Knowing the value of /3, one is in a position to calculate the transitional approximation to the 
overdense eigenvector 
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The Airy function approximation for the outer turning-point remains the same as before. Figure |l6| 
shows the first free eigenvector made up of the overdense barrier and Airy function approximations. 
20. Conclusion 

The diffraction of a plane wave by a sinusoidal potential is conveniently described by the Raman- 
Nath equation (Mathieu equation in Fourier space). A method for calculating the eigenvalues 
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(characteristic values) and continuizcd eigenvectors of this differential difference equation in the 
short wavelength limit is given. Working in Fourier space circumvents some of the difficulties 
associated with the infinite number of turning points inherent in a periodic potential. 

WKB-type solutions to the Raman-Nath equation serve as a starting point. Eigenvalues then 
follow from a simple Bohr-Sommcrfeld relation. Whilst the WKB solutions to the Raman-Nath 
equation still contain divergences, they reveal that the Raman-Nath equation can be interpreted as 
describing a wave in a double well potential, for which simple uniform approximations — solutions 
without singularities — exist in terms of the parabolic cylinder functions. There are three situations. 
Firstly, the 'bound' eigenstates lying below the separatrix in coordinate space lie above the central 
barrier in the double well in Fourier space and so a single uniform approximation in terms of Hcrmitc 
polynomials suffices. Secondly, at or just below the coordinate space separatrix the eigenstates lie at 
or just above the central barrier in the double well in Fourier space, causing two new turning points 
to appear. A complete eigenfunction with no singularities can be constructed by smoothly sewing 
together two transitional approximations — a parabolic cylinder function and an Airy function (the 
double well in Fourier space is symmetrical so only two transitional approximations are required 
rather than three). Thirdly, above the coordinate space separatrix and consequently below the 
central barrier in Fourier space, the 'free' eigenstates are also given by matching a parabolic cylinder 
function and an Airy function. A perscription is given for the modification of the Bohr-Sommcrfeld 
rule for the eigenvalues near the separatrix. 

When considering the diffraction of waves by a sinusoidal potential, knowing the cigenfunc- 
tions of the potential allows one to propagate the incident wave for any interaction distance and 
hence investigate dynamical diffraction phenomena which go beyond the phase-grating/Raman- 
Nath approximation, such as caustics (natural focussing). Scmiclassically, the superposition of 
eigenfunctions giving the diffraction pattern due to an incident plane wave only contains an ex- 
ponentially small contribution from the free eigenstates. However, as fi — > (A — > oo), even the 
number of 'bound' states become infinite. In a companion paper (O'Dell 2001) it will be demon- 
strated that a Poisson resummation of the eigenfunction superposition produces a new series each 
term of which is associated with classical paths belonging to a different topological class. Further- 
more, the number of terms required in this new sum depends linearly on the distance propagated 
through the potential (independent of the size of K) — only a finite number of terms are required 
for finite propagation distances, and so is computationally superior to the original eigenfunction 
sum which requires an infinite number of terms in the h —> limit. 
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Appendix: the method of uniform approximation 

Details can be found in the review by Berry and Mount (1972). The objective is to obtain an 
approximate solution of the Hclmholtz equation 



dq 2 



+ x(q)i>(q) = o 



in terms of solutions to one of the 'studied' equations, which will be written 

+ T{a)(f){a) = 0. 



da 2 



(99) 



(100) 



The choice of studied equation is determined by T(a) (not the Gamma function) being in some 
way similar to x{q)- This similarity implies that 4>(a) also resembles the wavefunction ijj{q) , and 
"can be changed into it by stretching or contracting it a little and changing the amplitude a little" . 
And so ip(q) "vvill be expressed in terms of <j){a) 



1>in) = f(q)H°(q))- 



Substitution of this definition into (fM) and making use of (100) leaves 



d 2 f 



dq 2 



da 



dq 



df da 



d 2 a 



da V dq dq 



dq 2 



0. 



(101) 



(102) 



The amplitude f(q) is as yet unspecified, so it is chosen to simplify (102) as much as possible. 
Putting 



/ 



renders ( |102[ ) into an equation purely for the 'mapping function' a{q) 



da 



dq 



x = x r- - — - 



da 



dq J dq 2 \ dq 



da 



(103) 



(104) 



which, when solved, gives a as a function of q. If a good choice of comparison function T(a ) ha s 
been made, then a(q) will be a slowly varying function and the second term on the rhs of (|l02| ) 
will be much smaller than the first. Clearly the criterion for this to be the case is 



1 





* d 2 




(1) 


dq 2 





<q)- 

When this is satisfied, the mapping relation reduces to 

x(q) 



< 1. 



da 
dq 



T(a) 



(105) 



(106) 



which through definition (103) also gives the amplitude /. Thus, by picking two points do and qo 

(107) 



which are 'equivalent', one finds a(q) from 

f yf±T(o~jda= f y/±ffij 
J an J an 



dq 



where the + or the — version can be chosen depending on the situation. The approximate solution 
to (|99|) is then 



x(q) 



(108) 
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In order for the comparison method to be viable, the mapping from q to a must be one to one, 



which requires that da/dq is never zero or infinite. Examining (106) this means that \ an d T must 
not diverge — which is assumed to be the case — and more relevantly, their zeros must be made to 
correspond. The zeros are of course the turning-points, and so, as Berry and Mount emphasise, 
"in the semiclassical limit all problems are equivalent which have the same classical turning-point 
structure''' . 

Perhaps the best known example of the uniform approximation is for the lone, first order (that 
is, the potential is locally linear) turning-point leading to the comparison equation 

0-^ = (109) 

whose solution is the Airy function, Ai(cr). Many potentials of interest are linear close to the 
turning point. As one moves away from the turning point the Airy function can be smoothly 
matched onto a WKB solution which is capable of handling very complicated potentials provided 
there are no turning points. Used in this way, as a patch across the turning point, the Airy function 
constitutes what is sometimes referred to as a transitional approximation. 
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A 


3 


fully numerical 


single well calc. 


modified calc. 


12500 


200 


0.996129 


0.996824 


0.996131 




198 


0.987197 


0.987337 


0.987199 




196 


0.976711 


0.976759 


0.976713 




194 


0.965430 


0.965461 


0.965432 




192 


0.953575 


0.953599 


0.953578 




190 


0.941244 


0.941264 


0.941247 




188 


0.928498 


0.928515 


0.928499 




186 


0.915379 


0.915394 


0.915381 


250000 


900 


0.999954 




0.999954 



Table 1: The bound eigenvalues near the separatrix: comparison of numerical result with 
standard Bohr-Sommerfeld condition for a well (|30|), and the modified conditions (|77j)-(f7q)- 
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Figure 1: A typical experimental set-up used in the investigation of atomic diffraction. 
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Figure 2: A series of plots showing the two momentum functions, p\ and p 2 , as functions of y for 
different values of (3. The top left has (3 = —0.999, each successive picture has (3 increasing by 0.2 
until the bottom right which has (3 = 1.201. It is the pi curve that dips down to zero when (3=1. 
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Figure 3: A comparison of the numerical (dots) th Bloch wave, out of 200 bound states, of the 
R-N matrix ( |To| ) for A = 12500, with its uniform approximation (solid line). (3 = —0.9937. Since 
the Bloch wave is symmetrical about y — 0, only the positive half is shown. 
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Figure 5: A comparison of the numerical (dots) 110 th Bloch wave, out of 200 bound states, of the 
R-N matrix ( |To| ) for A = 12500, with its uniform approximation (solid line). (3 = 0.2616. 
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Figure 7: A comparison of the numerical (dots) 200 th Bloch wave, out of 200 bound states, of the 
R-N matrix ( |To| ) for A = 12500, with its uniform approximation (solid line). (3 = 0.9961. 
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Figure 8: A comparison of the farfield wavefunction obtained by numerical diagonalisation 
(dashed), with the uniform calculation (solid), for A = 12500 and ( = tt/2. 
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Figure 9: A comparison of the farfield wavefunction obtained by numerical diagonalisation 
(dashed), with the uniform calculation (solid), for A = 12500 and ( = n. 
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Figure 10: A comparison of the farfield wavefunction obtained by numerical diagonalisation 
(dashed), with the uniform calculation (solid), for A = 12500 and ( = 3tt/2. 
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Figure 11: A comparison of the farficld wavefunction obtained by numerical diagonalisation 
(dashed), with the uniform calculation (solid), for A = 12500 and £ = 7n/2. 
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Figure 12: A comparison of the farfield wavefunction obtained by numerical diagonalisation 
(dashed), with the uniform calculation (solid), for A = 12500: ( = 81n/2. 
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Figure 13: A: The original p2, and B: transformed p2, phase momenta for f3 — 0.95. 




Figure 14: The WKB approximation using the undcrdense parabolic barrier action, see Equation 
(|7l|), and the fully numerical solution. In particular, this tests the derived phase angle /i as given 
by (f73|). The dots are the numerically calculated points, and the continuous line joins the WKB 
amplitudes. The value of A is 12500 and = 0.9961. 
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Figure 15: The parabolic barrier transitional approximation: A) as given by Equation @; B) 
the renormalised version. The dots are the fully numerical solution. The value of A is 12500 and 
(3 = 0.9961. 
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Figure 16: The 201 st Bloch wave, which is the first 'free' eigenvector. This requires both the 
overdense parabolic barrier solution, and an Airy function as transitional approximations. The 
two are joined at the 83 rd diffracted beam, which is at y — 0.742. The dots are the purely 
numerical calculation. 
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